program checksimul
use normeanvar
use datahadler
use matrixfuns
use imhoff
use linear_operators
implicit none
integer, parameter:: n=3,T=5000,num=100
integer, parameter:: p=n+(n*(n+1)/2)
double precision,parameter:: si=.999d0
double precision dmut(p,n),dvsig(p,n*n),sigma(n,n),matA(p,p),matB(n+1,p)&
&,matC(n+1,p+n+1),mtest1(n+1,1),vtest1(n+1,n+1),MD(p+n+1)&
&,VD(p+n+1,p+n+1),dmat(n,n),wmat(n+1,n+1),pr,dchiin,crit,b(n),poweru(num),powerc(num),eta(num)&
&,powerk(num),meank,vark,critk,dnorin,dnordf
integer i,j,k,ifault

eta=linspace(.00001d0,.02d0,num)
b=vassign(n,0.0d0)

forall(i=1:n+1,j=1:n+1)
	wmat(i,j)=0.0d0
end forall
forall(i=1:p,j=1:n)
dmut(i,j)=0.0d0
end forall
forall(i=1:p,j=1:n*n)
dvsig(i,j)=0.0d0
end forall
forall(i=1:n)
	dmut(i,i)=1.0d0
end forall
forall(i=1:n,j=1:n)
	dmat(i,j)=0.0d0
end forall
k=0
do j=1,n,1
	do i=1,n,1
		k=k+1
		dmat(j,i)=1.0d0
		dmat(i,j)=1.0d0
		dvsig(n+1:p,k)=vech(dmat)
		dmat(j,i)=0.0d0
		dmat(i,j)=0.0d0
	end do
end do
!sigma(1,:)=(/1.0d0,0.1d0,0.2d0/)
!sigma(2,:)=(/0.1d0,1.0d0,0.3d0/)
!sigma(3,:)=(/0.2d0,0.3d0,1.d0/)


sigma(1,:)=(/1.0d0,0.5d0,0.5d0/)
sigma(2,:)=(/0.5d0,1.0d0,0.5d0/)
sigma(3,:)=(/0.5d0,0.5d0,1.d0/)
!sigma=eye(3)
crit=dchiin(.95d0,dble(n+1))
critk=dnorin(.975d0)
do i=1,num,1
call meanvar(matA,matB,MD,VD,b,eta(i),si,sigma,dmut,dvsig,n,p)
matC(:,1:p)=matmul(matB,(.i.matA))
matC(:,p+1:p+n+1)=eye(n+1)

mtest1(:,1)=MD(p+1:p+n+1)*dsqrt(dble(T))
vtest1=(matC.x.VD).xt.matC

wmat(1,1)=2.0d0/(n*(n+2.0d0))
wmat(2:n+1,2:n+1)=(.5d0/(n+2.0d0))*(.i.sigma)


call imhofprob(pr,ifault,mtest1,vtest1,wmat,crit)
poweru(i)=1-pr
call imhofprob(pr,ifault,mtest1,VD(p+1:p+n+1,p+1:p+n+1),wmat,crit)
powerc(i)=1-pr
meank=mtest1(1,1)*dsqrt(wmat(1,1))
vark=vtest1(1,1)*wmat(1,1)
powerk(i)=1.0d0-dnordf((critk-meank)/dsqrt(vark))+dnordf((-critk-meank)/dsqrt(vark))
print*, i,ifault,poweru(i) !dot_product(mtest1(:,1),mtest1(:,1)),dot_product(mtest1(2:n+1,1),mtest1(2:n+1,1))
end do

call exportmat(eta,poweru,"p_aym3u.txt")
call exportmat(eta,powerc,"p_aym3c.txt")
call exportmat(eta,powerk,"p_aym3k.txt")

end program checksimul